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Biological systems are driven by intricate interactions among the complex array of molecules 
that comprise the cell. Many methods have been developed to reconstruct network models that 
attempt to capture those interactions. These methods often draw on large numbers of measured 
expression samples to tease out subtle signals and infer connections between genes (or gene prod¬ 
ucts). The result is an aggregate network model representing a single estimate for edge likelihoods. 

While informative, aggregate models fail to capture the heterogeneity that is often represented in a 
population. Here we propose a method to reverse engineer sample-specific networks from aggregate 
network models. We demonstrate the accuracy and applicability of our approach in several datasets, 
including simulated data, microarray expression data from synchronized yeast cells, and RNA-seq 
data collected from human subjects. We show that these sample-specific networks can be used to 
study the evolution of network topology across time and to characterize shifts in gene regulation 
that may not be apparent in the expression data. We believe the ability to generate sample-specific 
networks will revolutionize the field of network biology and has the potential to usher in an era of 
precision network medicine. 


1. INTRODUCTION 

In recent years we have come to appreciate that in 
many instances, and in particular in analyzing complex 
traits and diseases, a single gene or pathway cannot 
fully explain a particular cellular state. Instead, bi¬ 
ological processes are better characterized by complex 
networks whose structures are altered as the phenotype 
changes. Studying the pattern of connections between 
biological components, and how these structures change 
between cell states, can yield new insights into the mech¬ 
anisms driving disease. However, accurately reconstruct¬ 
ing these networks in a way that captures both the prop¬ 
erties and complexities of each phenotype remains a sig¬ 
nificant challenge. 

Estimating an informative network often relies upon 
combining information from large quantities of data 
(most commonly gene expression data) representing a 
potentially very diverse underlying sample population. 
This means that even in datasets where samples may 
represent a spectrum of phenotypes, most widely-used 
network inference methods use gene expression data to 
estimate an “aggregate” network Baca. These methods 
generally begin by calculating a statistic for each gene 
pair based on shared information across a set of input 
samples. These scores are frequently augmented in some 
way and are then used to infer the presence or absence 
of edges between the genes in the network. In this con¬ 
text, the heterogeneity in the underlying samples is gen¬ 
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erally seen as essential to correctly estimating a (single) 
network model given the variance in the data. However 
this consensus averaging may ignore the multiple differ¬ 
ent underlying regulatory networks that might be active 
in the individual samples. Recognizing this diversity is 
important. It is crucial that we examine the collection of 
networks to fully appreciate the diversity of processes at 
work in biological systems. 

Here we propose a method to reverse engineer sample- 
specific networks from an aggregate network model. The 
basic approach depends on comparing a network deduced 
from a population to a related network for the popu¬ 
lation missing a single sample and then estimating the 
contribution of that sample. What this means is that 
we can use aggregate network models to “extract” spe¬ 
cific networks for each of the individual input samples. 
We call this approach LIONESS (Linear Interpolation to 
Obtain Network Estimates for Single Samples). As will 
be shown in the examples below, this approach allows 
us to characterize the regulatory processes active in in¬ 
dividual samples and to model how regulatory networks 
change over time, both of which have many potential ap¬ 
plications in precision medicine and health and biomed¬ 
ical research. LIONESS is based on relatively straight¬ 
forward linear algebra and is independent of the network 
inference method used, suggesting the approach can be 
generalized to other single sample inference problems. 


1.1. LIONESS: Linear Interpolation to Obtain 
Network Estimates for Single Samples 

We begin by suggesting that an “aggregate” network 
predicted from a set of N samples can be thought of as 
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the average of individual component networks reflecting 
the contributions from each member in the input sample 
set. Mathematically, this means that the weight of an 
edge, between two nodes (i and j ) in an aggregate 
network derived using all samples (a) can be modelled as 
the linear combination of the weight of that edge across 
a set of networks, where each network in the set cor¬ 
responds to a sample used to reconstruct the aggregate 
model, and represents the relative contribution of 
sample (s) in that model: 

N N 
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Similarly, the edge weight in a second network recon¬ 
structed using all but one of the samples (a — q ), can be 
written as: 


N N 
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By comparing these two aggregate network models, we 
can solve for the network for an individual sample q. 
Subtracting the above equations we find: 
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The network specific to sample q in terms of the aggregate 
networks is then: 
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In summary, we find that the edge scores for a given in¬ 
dividual network (q) are equal to the difference in edge 
scores for an aggregate network constructed using all the 
samples, and an aggregate network reconstructed using 
all but the sample of interest, multiplied by a scaling fac¬ 
tor and added to the edge scores of the network recon¬ 
structed using all but the sample of interest (Figure [I]). 
In the following analysis we give samples equal weight 
(wq a ^ = 1 /N) although one could, in principle, weight 
samples differently based on the quality of the data for 
individual samples or some other measure. A more de¬ 
tailed version of the LIONESS derivation is provided in 
the Supplemental Material (Section [6|. 


2. RESULTS 

2.1. LIONESS accurately and reproducibly 
predicts networks using in silico data 

To test LIONESS we created a dataset where the un¬ 
derlying networks corresponding to each input expression 
sample are known (Supplemental Figure]]]). To begin, we 
created 10000 “gold standard” networks with 40 nodes 
and random edges but the same degree distribution. For 
each network, we generated 1000 random initial states 
and applied a Boolean model (see Methods) to deter¬ 
mine the corresponding attractors [3]. We then averaged 
over the all states defined within these 1000 attractors 
to generate “expression” values for the 40 nodes (which 
represent genes) in the network. 

We used this in silico data to evaluate the ability of 
LIONESS to accurately predict an individual sample net¬ 
work and explore how this prediction might be influ¬ 
enced by the number of “background” samples included 
when reconstructing the aggregate models. To begin, 
we used Pearson correlation to calculate aggregate net¬ 
works based on the in silico expression data. Each of 
these models is a complete graph with edge-weight esti¬ 
mates for every possible network interaction. We then 
applied LIONESS (Equation [29]) to reconstruct all indi¬ 
vidual sample networks, and evaluated each of the pre¬ 
dicted networks by comparing with the original “gold 
standard” networks and calculating the Area Under the 
Receiver Operator Characteristic curve (AUCROC, or 
more simply AUC). We repeated this multiple times and 
using a varying number of “background” samples. We 
observe that as we include larger numbers of samples in 
the aggregate network model, the accuracy of the single¬ 
sample networks predicted by LIONESS increases, and 
that these are consistently more similar to the “true” 
network than the aggregate models from which they were 
derived (Figure [2jA) . 

We also tested whether the LIONESS networks are re¬ 
producible. To do this we compared two single-sample 
networks that represent the same expression sample, but 
which were derived from two aggregate network mod¬ 
els that incorporate fully independent sets of background 
samples. We observe increasing reproducibility as we in¬ 
crease the number of background samples used to recon¬ 
struct the aggregate models, with near identity for large 
background sample sizes (Figure [2^3). As a control, we 
also compared two single-sample networks representing 
different gene expression samples that were derived us¬ 
ing the same set of background samples in the aggregate 
model. We see almost no similarity between these net¬ 
works, especially with increasing numbers of background 
samples. 

Finally, we investigated the generalizability of 
LIONESS by estimating single-sample networks from ag¬ 
gregate models derived using several common network 
reconstruction approaches, including Pearson correla¬ 
tion co-expression, PANDA [4], mutual information, and 
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Figure 1: Visual illustration of how LIONESS estimates the network for a single sample based on two aggregate network models, 
one reconstructed using all biological samples in a given dataset, and one using all except the sample of interest ( q ). 


CLR [5] (for more information, see Methods). We show 
the distribution in AUC values for the single-sample net¬ 
work predictions (Figure |2]E) and the difference in this 
AUC value compared to the AUC of the aggregate model 
evaluated using the same “gold standard” (Figure [2p). 
We also report the median of these values and the percent 
of single-sample networks whose AUC is greater than the 
AUC of the overall aggregate model (Figure [2]G— D). We 
find that LIONESS is able to accurately predict single¬ 
sample networks for all four of the approaches, with me¬ 
dian AUC values between 0.6 and 0.7 (Figure [2p). In 
all cases the majority of the single-sample networks re¬ 
constructed using LIONESS are more predictive of the 
true underlying networks than the aggregate from which 
they were derived. Of the four network inference meth¬ 
ods, PANDA networks were the most specific, with over 
99% of the single-sample network models more represen¬ 
tative of the true underlying network than the aggregate 
network (Figure [2jD). 


2.2. Estimating single-sample networks using 
experimental data from yeast 

We next tested LIONESS using experimental data 
from cell-cycle synchronized yeast cells. We downloaded 


gene expression data (GEO accession, GSE4987; [6]) con¬ 
sisting of dye-swap technical replicates measured every 
five minutes for 120 minutes and ma-normalized each 
replicate [7:; we then removed probe sets with missing 
data, batch-corrected using ComBat 0, averaged probe 
sets mapping to the same ORF annotation and quantile- 
normalized the resulting gene-by-sample matrix of ex¬ 
pression values. We note that the 105 minute time point 
was excluded in both replicates due to poor hybridization 
performance [6]. 

We used four different approaches (Pearson Correla¬ 
tion, PANDA [4], mutual information, and CLR 0 ) to 
reconstruct aggregate networks for this dataset and ap¬ 
plied LIONESS to estimate the networks for each of the 
individual samples. The correlation between each pair 
of the estimated sample-specific networks is shown in the 
first column of Figure[3^ (Rl&R2-from-Rl&R2). We see 
that network estimates for the same technical replicate 
are highly similar, as evidenced by the strong diagonal 
in the upper-right and lower-left square of each compar¬ 
ison; additional structure is also evident in off-diagonal 
similarities that reflect the fact that the time course data 
represents more than one cell cycle. 

To test if strong reproducibility was due to including 
replicates in the expression data, we also ran LIONESS 
separately on each individual replicate. This analysis 
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Figure 2: Evaluation of LIONESS’ ability to recover known single-samples networks in an in silico dataset. (A) A plot of the 
mean and standard deviation of AUC values calculated when estimating single-sample networks from aggregate models that 
include information from varying numbers of background samples. (B) The median and inter-quartile range of the similarity 
(calculated using the Spearman correlation) between pairs of networks, as a function of the number of background samples 
used to calculate the aggregate model. For the plots in (A) and (B) each point represents the results across 1000 random 
selections of samples. (C) The distribution of AUC values that result when we evaluate each of the 10000 LIONESS-predicted 
single-sample networks using their corresponding original “gold-standard” networks. The values are almost always greater 
than 0.5, indicating overall accurate prediction of the true sample-specific networks. (D) The difference in AUC values when 
evaluating the 10000 LIONESS-predicted networks versus evaluating the aggregate network reconstructed using all samples. 
This difference is almost always greater than zero, indicating that the LIONESS-predicted single-sample networks are more 
reflective of the underlying true networks than the aggregate network from which they are derived. 


produced 24 single-sample networks estimated using only 
the data in replicate one, and 24 single-sample networks 
estimated using only the data in replicate two (Rl-from- 
R1 & R2-from-R2). The correlation between these net¬ 
works is shown in the second column of Figure [3jA. As 
before, we observe strong reproducibility between tech¬ 
nical replicates. Interestingly, even though we have cor¬ 
rected for batch effects in our expression data, several of 
the methods, especially CLR, appear to be sensitive to 
the “background” data used. 

Although we have shown the networks estimated by 
applying LIONESS to each of four reconstruction ap¬ 
proaches are reproducible , it is also important that these 
networks are also accurate. To test the accuracy of each 
of the four network reconstruction methods, we compared 
the aggregate networks to ChIP-chip information, and 
tested whether the aggregate network models could pre¬ 
dict transcription factor binding. For three of the four 
network reconstruction approaches we observe AUC val¬ 
ues close to 0.5 (Figure |3j3). The exception is PANDA, 
which had an AUC above 0.7. 

PANDA differs from the other approaches in that it 
begins with a “prior” gene regulatory structure (based 
on motif data; see Methods) and uses the Pearson cor¬ 
relation as a proxy to gene co-regulation by the same 
upstream regulator, rather than a direct regulatory in¬ 
teraction between a regulator and a downstream target. 


We also note that, in addition to having overall higher 
AUC values, the improvement in AUC values compared 
to a null state (the AUC of the original “seed” motif prior 
for PANDA and AUC= 0.5 for the other approaches) is 
greatest for PANDA (Figure [3^). Based on these results, 
in the following analysis we focus on the single-sample 
networks derived using PANDA as the underlying ag¬ 
gregate approach. Results for the other reconstruction 
approaches are presented in Supplemental Figure [2] 


2.3. Single-sample networks show periodic 
structure across the cell cycle 

Next we averaged sample networks representing the 
same time point in each of the two replicate network re¬ 
constructions and compared the results to the ChIP-chip 
data to evaluate the accuracy of each of the resulting 
network models. We recognize that this is not a true 
evaluation of the sample-specific networks, since we ex¬ 
pect the network structure to change throughout the cell 
cycle. In Figure [4]A we show the difference in the AUC 
values for the each of the PANDA+LIONESS-predicted 
single-sample networks when compared to the AUC of 
the motif prior network that was given as an input to 
PANDA. Although the increases are modest, each of the 
single-sample networks has an overall AUC greater than 
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Figure 3: Analysis of LIONESS networks predicted for 48 ex¬ 
pression samples collected across a yeast cell-cycle time course 
experiment. (A) LIONESS was used to predict networks for 
each sample in the expression dataset by applying four differ¬ 
ent aggregate network reconstruction approaches. For each 
approach we built the aggregate models either using all sam¬ 
ples (R1&R2 from R1&R2), or only the samples from the 
same technical replicate (Rl-from-Rl & R2-from-R2). The 
Spearman correlation was used to evaluate how similar these 
networks are to each other. (B) A bar graph showing the 
change in AUC values for the aggregate models predicted for 
each approach, compared to a baseline. For Rl&R2-from- 
R1&R2 (red bars) we have a single aggregate network. For 
Rl-from-R2 &; R2-from-R2 we have two aggregate models 
(one for each replicate), which we average before comput¬ 
ing the AUC (blue bars). In this analysis, the baseline AUC 
values for the Pearson correlation, MI and CLR networks was 
set equal to 0.5. Since PANDA begins with a prior network 
derived from motif structure, we use that AUC value (0.68) 
as the baseline in order to have a more fair comparison with 
the other approaches. 


the prior network. There appears to be an oscillatory 
pattern in these AUC values across the time course, with 
higher values at greater times post the initial synchro¬ 
nization. This makes sense as cells will gradually become 
unsynchronized, better representing the cell population 
used for the ChIP-chip data. 

We next tested whether these single-sample networks 
could provide insight into gene regulation and dynamic 
cellular network properties. We identified the 1000 edges 
with the highest variability across the individual net¬ 
works and plotted those edge weights in Figure [4^. We 
again observe strong oscillatory patterns that correspond 
closely with the oscillations in AUC. It is worth noting 


that all these highly variable edges originate from one of 
four transcription factors, each of which is important in 
regulating the yeast cell cycle (MBP1, SWI4, SWI6, and 
STB1) 0. 

We examined the genes for which there is strong ev¬ 
idence of targeting by these transcription factors (aver¬ 
age edge weight across all LIONESS networks greater 
than zero). In Figure [4]C we plot the average weight of 
these high-evidence interactions for each regulating tran¬ 
scription factor and the average expression of their tar¬ 
get genes. It is immediately apparent that oscillation in 
edge weights occurs at exactly twice the frequency of the 
oscillation in gene expression, and that the expression os¬ 
cillation has a period approximately equal to the yeast 
cell cycle. 

To understand this result we point out that PANDA 
interprets correlation in target gene behavior as indica¬ 
tive of upstream regulation, emphasizing edges between 
a regulator and its targets either when those targets are 
all coordinately expressed (activated), or coordinately 
not-expressed (de-activated or repressed). Thus high 
edge weight can be interpreted as evidence for informa¬ 
tion flow, which may be the result of the physical pres¬ 
ence of a transcription factor (TF) regulating its down¬ 
stream targets, but could also reflect a strong lack of that 
TF. Indeed, by using correlation-based metrics to infer 
network relationships many reconstruction approaches 
operate under a similar mathematical framework. Al¬ 
though this observation may change how we interpret the 
edges in the estimated single-sample networks, the “turn 
on/turn off” behavior is exactly what one would predict 
given how PANDA estimates networks and is further evi¬ 
dence that LIONESS is extracting realistic single-sample 
networks. 


2.4. Reconstructing single-sample networks for 
human lymphoblastoid cell-lines 

Lastly, we wanted to test LIONESS in the context of 
modeling human gene regulatory networks. Based on 
the results we observed when applying four different ag¬ 
gregate reconstruction approaches in conjunction with 
LIONESS (see Figures[2jj3]), we chose to focus on PANDA 
as our aggregate network reconstruction approach in this 
analysis. PANDA incorporates prior network informa¬ 
tion, which we constructed by using Haystack [10] to scan 
the promoter regions of human genes for JASPAR tran¬ 
scription factor binding motifs m ■ 

For our expression data we used a set of RNA- 
seq experiments performed across multiple immortal¬ 
ized lymphoblastoid cell lines representing different in¬ 
dividuals [12]. We downloaded raw fastq files from the 
Pritchard lab website (http://eqtl.uchicago.edu/) 
and aligned samples to hgl9 using Bowtie m-, subse¬ 
quent quality control analysis using RNA-SeQC [14] ex¬ 
cluded two individuals from further analysis due to low 
expression profile efficiency scores. This left us with a 
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Figure 4: Characterizing networks across the yeast cell-cycle. (A) The difference in the AUC for each sample-specific network 
predicted by PANDA+LIONESS versus the AUC of the motif prior used to seed the PANDA reconstruction approach. The 
same ChIP-chip data was used as the basis for each evaluation. (B) A heat map of the edge weights for the 1000 most variable 
edges across the sample-specific network models. Rows are Z-score normalized for visualization purposes. (C) The average 
expression of genes targeted by the four transcription factors that were identified as regulatory nodes of the 1000 top most 
variable edges as well as the average weight of high-confidence edges that extend between those transcription factors and their 
target genes. 


final set of 153 RNA-seq experiments that includes repli¬ 
cates and represents 65 distinct individuals. We normal¬ 
ized this data using DEseq2 m- For additional data pro¬ 
cessing and normalization information, see Online Meth¬ 
ods. 

We used PANDA to reconstruct aggregate gene regu¬ 
latory network models that integrate transcription factor 
binding information and expression data. We then ap¬ 
plied LIONESS to these aggregate models, resulting in 
153 single-sample networks, one for each of the RNA- 
seq expression samples. A hierarchical clustering (com¬ 
plete linkage, Spearman Correlation) of the network edge 
weights demonstrates that networks for the same indi¬ 
vidual nearly always cluster more strongly with each 
other than with networks representing different individ¬ 


uals (Supplemental Figure [3]). We also evaluated these 
networks using a ChIP-seq “gold standard” and found 
91.5% had a higher AUC compared to the “seed” motif 
network used by PANDA (see Online Methods). This 
analysis demonstrates that even when constructing net¬ 
works using biological data from higher-order organisms 
such as human, the sample-specific networks predicted 
by LIONESS are highly reproducible and predictive of 
known interactions. 
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2.5. Complex relationships between network 
targeting and gene expression 

In our previous analysis in yeast (Figure J4p) we ob¬ 
served that edges in single-sample networks reflect strong 
information flow, which does not necessarily correspond 
to physical binding events. To explore this further, we 
investigated the relationship between gene targeting and 
expression in human networks. First, we averaged single¬ 
sample networks that represent the same individual, re¬ 
sulting in 65 “person-specific” regulatory networks. We 
then selected high-evidence regulatory interactions for 
each transcription factor (average edge-weight across all 
single-sample networks greater than zero), and directly 
compared the mean edge-weight for these interactions in 
each of the single-sample networks to the average ex¬ 
pression of the targeted genes in the original expression 
samples. 

We consistently observe a non-linear relationship be¬ 
tween targeting and expression, with the highest aver¬ 
age edge weights occurring when target genes have ei¬ 
ther overall high or low expression levels (Figure I K); 
this is consistent with what we observed in our yeast 
analysis (Figure |4p). Coloring by the transcription factor 
expression level in each sample reveals further patterns; 
some transcription factors appear to be acting primarily 
as activators (increased target gene expression upon in¬ 
creased TF expression and targeting) or as repressors (de¬ 
creased target gene expression upon increased TF expres¬ 
sion and targeting). However, the relationship between 
transcription factor and target genes is not always sim¬ 
ple, indicating that other regulatory mechanisms, such as 
co-activators, post-translational modifiers, or epigenetic 
mechanisms, are likely playing an important role in me¬ 
diating these regulatory events. 


2.6. Increased network targeting corresponds to 
open chromatin 

There is available DNase hypersensitivity profiling of 
cell lines representing these same 65 individuals [Jjj. As 
we did not use this data when modeling the LIONESS 
networks, we next sought to leverage this information to 
investigate how network structures might reflect epige¬ 
netic state. To begin, we downloaded the data from the 
Pritchard lab website (http://eqtl.uchicago.edu/) 
and called DNase “peaks” for each sample. When a peak 
fell within the promoter region of a gene, we associated 
that gene-sample pair with a score reflecting the signifi¬ 
cance level of the associated peak call. We found 12424 
genes with a DNase promoter-peak in at least one sample 
and 3488 with a promoter-peak in all samples. For more 
information on the DNase data processing, see Online 
Methods. 

A DNase hypersensitivity peak represents a region of 
open chromatin that is often presumed to be occupied by 
a set of regulatory proteins, including transcription fac¬ 


tors. Consequently, we wanted to determine if differences 
in chromatin state between the 65 individuals is reflected 
by alterations in transcription factor targeting across our 
single-sample networks. Based on the analysis presented 
in Figure [4p and[5]4, we hypothesized that physical bind¬ 
ing of a transcription factor to the promoter region of a 
target gene would be best represented by combining: (1) 
the edge weight of that interaction in our single-sample 
network models, since this value indicates whether infor¬ 
mation is flowing between that transcription factor and 
target gene; and (2) the expression level of the transcrip¬ 
tion factor, since this value indicates whether the TF is 
physically present in the cell (Figure^). Therefore, we 
combined this information (see Online methods), result¬ 
ing in a set of expression-modified edge-weights for each 
sample. 

We used the sum of the edge-weights associated with 
each gene to estimate the number of transcription fac¬ 
tors regulating that gene in each of the 65 single-sample 
networks ( k ^ L+e ^). As a control we also calculated gene¬ 
targeting two other ways: (1) using LIONESS edge- 
weight estimates in the absence of gene expression in¬ 
formation ( k ( L )) and (2) using gene expression informa¬ 
tion in the absence of LIONESS-predicted edge-weights 
(fc( m+e )); for the second measure we combined transcrip¬ 
tion factor expression in each sample with the motif in¬ 
formation used for PANDA’s prior. Finally, to evalu¬ 
ate the association of network targeting with chromatin 
state, for each gene we calculated the Spearman correla¬ 
tion between these targeting-values across the networks 
and the significance scores of that gene’s promoter-DNase 
across the corresponding cell-lines. We find that gene¬ 
targeting in the expression-modified LIONESS model 
(^(L+e)) ^ ver y strongly correlated with promoter-DNase 
events (Figure[5]C). This association is greater than when 
using only expression and motif information (&( m + e )) 5 
demonstrating that the LIONESS approach provides ad¬ 
ditional information on chromatin state not apparent in 
the data used to seed the algorithm. 


2.7. Differential-targeting of genes highlights 
important biological processes 

Finally, we wanted to determine if there are common 
structures across these single-sample regulatory networks 
that might be reflective of important biological processes. 
We performed a hierarchical clustering (complete link¬ 
age; Spearman Correlation) on the edge-weights in the 
65 single-sample networks and observe clear sets of sam¬ 
ples defined by this network information (Figure [6pA) . To 
evaluate if this clustering would also be identified using 
the expression data alone, we also performed a hierarchi¬ 
cal clustering using gene expression values (Figure [6^). 
We observe a strikingly different grouping of individuals 
compared to the one obtained when clustering the net¬ 
works. 

To determine if the expression-based clustering high- 
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Figure 5: Comparison of gene regulation, gene expression and DNase hypersensitivity data. (A) For six representative transcrip¬ 
tion factors, the mean expression of target genes and the mean weight of the edges targeting those genes across the 65 samples 
are plotted. For each sample, the expression of the TF is shown as a color, scaled to the normal distribution for visualization 
purposes. (B) A cartoon illustrating how high edge weights and thus regulatory activity is not necessarily equivalent to the 
presence of a physical interaction. (C) The distribution of the Spearman correlation values when comparing gene-targeting 
(calculated by combining LIONESS predictions with TF expression; U L+e ^) and the significance level of DNase hypersensi¬ 
tivity in a gene’s promoter across all the samples. We also show the percentage of genes whose targeting positively correlates 
with DNase hypersensitivity when targeting is calculated using only the LIONESS-predicted edge weight ( k no expression 
considered) or a combination of expression and motif information ( k ( m + e )) > These percentages were calculated either using 
all 12424 genes included in our network model, or for the set of 3488 genes with a DNase-peak called in all 65 samples (open 
genes). 


lights important biological information we compared the 
expression levels of genes between groups of individu¬ 
als and performed Gene Set Enrichment Analysis m • 
Interestingly, although we find many genes that are dif¬ 
ferentially expressed (2620 with FDR < 0.01), running 
GSEA finds no enrichment for known biological func¬ 
tions. Separately, we ran GSEA to evaluate if biolog¬ 
ical functions are associated with differential-targeting 
patterns between the two groups of networks shown in 
Figure m • This analysis identified enrichment for 
many cellular processes related to cell proliferation (in 
the smaller “orange” cluster; n = 18) and immune func¬ 
tion (in the larger “blue” cluster; Figure J6p). 

Unfortunately, there is little publicly available pheno¬ 
typic information for the 65 individuals in our dataset, 
and the characteristics that we do know m are not 


highly associated with either the expression or the net¬ 
work clusters. However, given our functional enrichment 
results, we believe that the regulatory differences we ob¬ 
serve between the network clusters is likely related to 
differences in cellular growth rate induced by variable 
Epstein-Barr Virus (EBV) levels in the cell-lines. EBV is 
used to transform human B-cells into immortalized lym- 
phoblastoids and is known to activate NF/H3 transcrip¬ 
tional response [15]. Consistent with this hypothesis, we 
find the signature “Activation of NF-kappaB in B-cells” 
highly targeted in the small, “cell proliferation” cluster 
(ES= 0.5, FDR= 2 * 10 -3 ). 

These results indicate that evaluating single-sample 
networks can lend insight into alterations in the bi¬ 
ological processes active in different individuals even 
when a similar analysis of the gene expression data does 
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Figure 6: (A) A hierarchical clustering on the edge-weights for 65 regulatory networks, one for each subject included in the 
RNA-seq dataset. (B) An equivalent hierarchical clustering on gene expression values for these 65 individuals. This clustering 
is distinct from the one based on the network edge weights. Subject labels are colored based on the network clustering. (C) 
Reactome pathways enriched based on a GSEA analysis using gene-targeting instead of gene expression and comparing the 
right and left clusters of networks presented in (A). No Reactome pathways were identified when comparing the expression 
values of genes in the different groups defined by the hierarchical clustering presented in (B). 


not. Further, the relevance of the processes identi¬ 
fied through differential targeting in the model suggests 
that LIONESS+PANDA are capturing essential biologi¬ 
cal processes relevant to each sample. 


3. DISCUSSION 

Network models that use aggregate statistics to de¬ 
scribe regulation across a population are important and 
have contributed to our understanding of a wide range 
of phenotypes. However they are limited to modeling 
the biological processes common across samples, and thus 
cannot fully characterize any potential heterogeneity in 
regulatory processes within set of underlying samples. 
This heterogeneity is important. If we think about the 
simplest system—a collection of cells from a single organ, 
we expect each individual cell to have its own unique gene 
expression levels and its own unique active gene regula¬ 


tory networks. Understanding the diversity and complex¬ 
ity of these networks is essential if we are to fully under¬ 
stand the population dynamics of the system. Similarly, 
in analyzing disease phenotypes, we recognize that each 
individual has his or her own unique manifestation of the 
disease and that this should be reflected in individual- 
specific gene regulatory networks. Our current aggregate 
network inference methods do not allow us to fully ex¬ 
plore the complexity of the systems we are hoping to 
understand. 

We propose LIONESS as a theoretical framework to re¬ 
construct networks at the level of the individual sample 
by taking advantage of the aggregate models predicted 
by widely-used network reconstruction methods. The 
LIONESS approach not only addresses the problem of 
phenotype heterogeneity, but also provides a means of 
estimating networks when samples of a particular pheno¬ 
type or disease subtype are rare. One can pool a diverse 
collection of samples and then extract individual sam- 
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pie networks, grouping these afterward based on network 
edge weights. 

One implicit assumption in the LIONESS frame¬ 
work is that the samples in a given collection repre¬ 
sent the full span of possible regulatory interactions. 
Namely, LIONESS assumes that an aggregate model re¬ 
constructed without the sample of interest represents a 
common network, and that the difference in this model 
with the one constructed using all samples represents a 
perturbation to that common network. Importantly, in 
our analysis we observe that increasing the number of 
samples used as a “background” when reconstructing an 
aggregate model is beneficial to the single-sample net¬ 
works. 

We recognize that the quality and interpretation of the 
LIONESS networks is highly dependent on the underly¬ 
ing aggregate network reconstruction algorithm used. We 
emphasize that there are many existing network recon¬ 
struction methods and there remains much debate in the 
field as to the “correct” or “best” one to use. In this 
paper we show only four representative gene expression 
network reconstruction approaches: Pearson correlation, 
Mutual Information, Context Likelihood of Relatedness, 
and PANDA. These were chosen not because they are 
necessarily the best, but because they illustrate network 
reconstruction methods that use either a linear (Pearson) 
or non-linear (mutual information) correlation measure, 
and the extensions of those measures to better capture 
true regulatory interactions instead of simple correlative 


effects. However, within this limited sample, our analysis 
suggests that applying LIONESS to aggregate networks 
reconstructed using PANDA has the greatest potential 
for reconstructing accurate network models able to iden¬ 
tify relevant network biology. Finally, although we tested 
our approach in the context of using gene expression to 
reverse-engineer regulatory networks, the linear frame¬ 
work at the heart of LIONESS is generalizable and could 
be applied to infer sample-specific networks of other types 
or using other data sources and network inference meth¬ 
ods. 

Looking ahead, we believe that accurately estimating 
single-sample networks provides a unique opportunity to 
investigate and identify patient-specific disease processes 
and has the potential to open the way for more person¬ 
alized network medicine. In network medicine applica¬ 
tions regulatory models are often constructed for a dis¬ 
ease subtype. Applying LIONESS to these models could 
identify the sub-networks that are active in an individual 
patient, providing an opportunity to select specific ther¬ 
apeutic interventions based on the unique characteristics 
of the gene regulatory networks active in each individual. 
Finally and critically for medical applications, we point 
out that by reconstructing patient-specific networks we 
will be able to directly correlate network responses and 
features with clinical information, opening the door for 
investigating the complexity of disease features within 
the context of network structures. 
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4. ONLINE METHODS 
4.1. Network reconstruction approaches 

Pearson Correlation: Pearson correlation evaluates the 
degree of a linear relationship between two variables. 
Regulatory networks can be reconstructed by calculat¬ 
ing the Pearson correlation coefficient between expression 
levels of each TF and each target gene. These coefficients 
are measures of whether TFs and target genes are being 
co-expressed, which may indicate a regulatory event. 

Passing Attributes between Networks for Data Assimi¬ 
lation (PANDA): PANDA builds regulatory networks by 
starting with a prior of possible interactions between TFs 
and target genes, for example TF motif binding infor¬ 
mation. PANDA integrates this regulatory prior with 
gene expression information and protein-protein interac¬ 
tion data, using a message passing approach to deter¬ 
mine information flow between the different data types. 
The message passing algorithm that is used in PANDA 
is based on the fact that if expression levels of two genes 
correlate, those genes are more likely to be regulated by 
similar sets of TFs than two genes that do not show corre¬ 
lation in expression. Similarly, TFs that can form com¬ 
plexes are more likely to regulate similar sets of target 
genes than TFs that cannot bind together. 

Mutual Information (MI): MI compares the joint prob¬ 
ability distribution of two variables to the products of 
their corresponding marginal distributions. Similar to 
Pearson correlation, MI is a measure of association be¬ 
tween TFs and target genes, which may indicate a regu¬ 
latory event. This method does not assume linearity or 
continuity of the data used for building the network. 

Context Likelihood of Relatedness (CLR): CLR is 
based on MI, but applies a double z-score transformation 
to the MI scores, normalizing each TF-gene interaction to 
the background distribution of MI values for each gene, 
and to the background distribution of MI values for each 
TF. 


4.2. Generation of the in silico expression data and 
regulatory networks 

To test LIONESS’s ability to reconstruct sample- 
specific network models we generated a set of network 
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models and a corresponding associated set of gene ex¬ 
pression profiles. To begin, we created a single “seed” 
network model with forty nodes. To approximate the 
structure of biological networks [20], the out-degree of 
nodes in this seed model were given a power-law degree 
distribution (generated using the approach published in 
Clauset et al. [21], with a = 5) with their targets selected 
randomly. We ensured that the out-degree and in-degree 
of all nodes in this seed network was greater than zero 
and less than forty. The resulting “seed” network had 81 
total edges and a density around 5%. 

Next we randomized this “seed” network model, hold¬ 
ing the degree fixed. Then we generated a set of initial 
Boolean states for each node in the network, and deter¬ 
mined the subsequent states of the nodes using Stouffer’s 
Z-score method: 




=round 


CDF- 1 



i 


( 5 ) 


where CDF 1 is the inverse cumulative distribution 
function for the normal distribution, Zij is the z-scored 

weight of an edge from node i to node j, and is the 
state of node i at time t. This Boolean model was run un¬ 
til an attractor solution was found. In total we generated 
1000 random initial Boolean states for each randomized 
network, resulting in 1000 attractor solutions. The ex¬ 
pression level of a node in the randomized network was 
then estimated as the average across these steady-state 
solutions. This entire process was then repeated for other 
randomized versions of the “seed” network model. In 
total we performed 10000 network randomizations and 
generated 10000 corresponding expression states. 


4.3. Processing the yeast cell cycle expression data 


GPR files were downloaded from the Gene Ex¬ 
pression Omnibus (GEO; accession GSE4987). Each 
of two replicates were separately ma-normalized us¬ 
ing the maNorm() function in the “marray” library in 
R/Bioconductor [7 . The data was batch-corrected us¬ 
ing the ComBat() function in the “sva” library [22] and 
probe-sets mapping to the same gene were averaged, re¬ 
sulting in expression values for 5088 genes across fifty 
conditions. Two samples (corresponding to the 105 
minute time point) were excluded for data-quality rea¬ 
sons, as noted in the original publication, and genes with¬ 
out motif information (see below) were then removed, 
giving a final expression data-set containing 48 samples 
and 3551 genes. This data was quantile-normalized and 
used in all subsequent analysis. 


4.4. Generating the yeast motif prior data for 
PANDA and yeast ChIP-chip gold standard 

PANDA requires a prior regulatory network structure 
in addition to gene expression information. To construct 
a motif prior network for yeast we downloaded predicted 
binding sites for 204 yeast transcription factors [23, [24]]. 
This data includes 4360 genes with tandem promoters. 
3551 of these genes are also covered on the gene expres¬ 
sion array (see above). 105 total transcription factors 
in this dataset target the promoter of one of these 3551 
genes. The motif map between these 105 transcription 
factors and 3551 target genes was used as a prior regula¬ 
tory network input to the PANDA algorithm. A subset of 
65 of these transcription factors that also had expression 
information was used to reduce the size of the Pearson, 
MI and CLR predicted networks to edges that extend 
between transcription factors and genes. 

We used ChIP-chip data from Harbison et al. [23] 
to calculate AUC values for the predicted single-sample 
yeast regulatory networks. Targets of transcription fac¬ 
tors in this dataset were defined using the criterion 

p < 0.001. 


4.5. Processing the human DNase hypersensitivity 

data 

Raw DNAse hypersensitivity data was down¬ 
loaded from the Pritchard lab website (http: 
//eqtl.uchicago.edu/dsQTL_data/RAW_DATA/, ac¬ 

cessed June 2014). 204 different samples corresponding 
to 70 different cell lines were available for download. 
Data were aligned to the hgl9/GRCh37 reference 
genome using Bowtie m, with options -v 1 and -m 
10, allowing for not more than one mismatch, and 
suppressing any alignments for reads having more than 
ten reportable alignments. Quality control using Bowtie 
output revealed that 16 samples had a high percentage 
(greater than 80%) of failed reads. We called DNase 
hypersensitivity peaks using MACS [25]. Peaks with 
significance score of less than 10 -5 were mapped to the 
nearest gene using HOMER [2^. When the peak fell 
within the promoter region of the gene, we assigned a 
score to that sample-gene pair equal to —10 * logio(p ); 
otherwise the sample-gene pair was given a default score 
of zero. We then removed the 16 poor quality samples, 
as well as samples for which we did not have good 
quality RNA-seq data (see below), leaving us with 177 
samples corresponding to 65 cell-lines. 

To obtain a promoter-DNase score specific for each cell¬ 
line, we averaged technical replicates. We then filtered 
this data to include only genes for which we also had 
RNA-seq data available (see below) and genes that were 
present in our motif prior. The result was a matrix of 
DNase-promoter scores that included 12424 genes that 
had a DNase promoter-peak in at least one sample (at 
least one row-entry greater than zero) and 3488 had a 
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promoter-peak in all samples (all row-entries greater than 
zero). 


4.6. Processing the human RNA-seq data 

RNA sequencing (RNA-Seq) data were down¬ 
loaded from the Pritchard lab website (http://eqtl. 
uchicago.edu/, accessed April 2014). 173 different 

samples corresponding to 74 different cell lines were 
available for download. We aligned all samples to the 
hgl9/GRCh37 reference genome using Bowtie [13], with 
options -n 3 and -m 1, allowing for not more than three 
mismatches in the seed (28 bases on the high-quality end 
of the read, the default in Bowtie), and suppressing all 
non-unique alignments. We used RNA-SeQC El to de¬ 
termine the quality of the reads, using an expression 
profiling efficiency cut-off of 0.75. Cell lines NA19119 
and NA18853 fell below this cut-off. Next we used sub¬ 
read m to count reads, and the subread algorithm fea- 
tureCounts to assign and summarize counts to genes. Fi¬ 
nally, we removed samples with poor quality reads, and 
samples for which we did not have good quality DNase 
hypersensitivity data available (see above), leaving us 
with 153 samples corresponding to 65 different cell lines. 

We used the DESeq2 [15] package to analyze read 
counts for these 153 samples and adjusted for different li¬ 
brary sizes using the estimateSizeFactors function. Only 
genes that had raw counts in at least 50% of all 153 sam¬ 
ples (21516/57820, or 37% of all genes) were retained for 
further analysis. Correction for gene length was also per¬ 
formed; each intensity value was divided by the length of 
the corresponding gene (defined as the total length of 
the genomic region covered by the features/exons) in the 
“gene_id” met a-feature in featureCounts. Finally, En- 
sembl gene ids were converted to HGNC gene symbols 
(16901/21516 genes) using R package biomaRt p8l [29] . 
This gene list was subsetted to only include genes for 
which we found at least one transcription factor binding 
site in the adjacent promoter (see below), and for which 
we had at least one DNase hypersensitivity peak-call in 
the adjacent promoter (see above). This resulted in a 
matrix of 12424 HGNC symbols by 153 expression sam¬ 
ples. 


4.7. Generating the human motif prior data for 
PANDA and human ChIP-seq gold standard 

We downloaded JASPAR motifs (http://jaspar. 
genereg.net/; El) and then used Haystack m to scan 
the entire hgl9 genome for these motifs. Of the 205 mo¬ 
tifs in the JASPAR database, only 158 had genomic hits 
that met our significance threshold (p < 10 -6 ). We used 
HOMER (http: //homer. salk. edu/homer/ngs/index. 
html; [26]) to get the distances of these motif hits to the 
nearest transcriptional start site (TSS). We used these 
reported distances to parse the motif hits based on their 


TSS proximity, keeping only those hits within the “pro¬ 
moter” region of a gene, which we define as [—750, +250] 
around the TSS. We then filtered for genes for which 
we had DNase hypersensitivity data and RNA-seq data 
available. This information was used to make a prior 
transcription factor to gene map that we used when run¬ 
ning PANDA. 

To generate a gold-standard to evaluate the predicted 
network models, we downloaded ChIP-seq bed files as¬ 
sessing the localization of 39 transcription factors which 
were included in our motif network and have been as¬ 
sayed in GM12878 by ENCODE (http: //genome. ucsc. 
edu/ENCODE/). We used this location information to cre¬ 
ate a “gold standard” transcription factor to gene net¬ 
work, using HOMER (http://homer.salk.edu/homer/ 
ngs/index .html; [26]) to map TF peaks to the nearest 
gene in a process identical to our mapping of TF-motif 
location information, described above. 


4.8. Running LIONESS on networks reconstructed 
using mutual information and CLR 

In order to calculate the mutual information and re¬ 
construct networks based on the Context Likelihood of 
Relatedness, we used the build.mim() and clr() func¬ 
tions (default parameters) within the “minet” package 
in R/Bioconductor [30] . It is worth noting that these 
algorithms create symmetric gene by gene matrices of 
predicted edge scores. Therefore, in the context of re¬ 
constructing single-sample yeast regulatory networks, we 
reduced these aggregate networks by selecting only the 
portion of the predicted gene-by-gene matrix that corre¬ 
sponds to edges from a transcription factor to a target 
gene. 


4.9. Running LIONESS on networks reconstructed 
using PANDA 

PANDA requires a prior regulatory network structure 
in addition to gene expression information. For the 
in silico data we used an identity matrix (correspond¬ 
ing to each transcription factor targeting only itself) as 
the prior regulatory network. For the yeast and human 
data a prior regulatory network was constructed based 
on transcription sequence-motif information (see above). 
It is worth noting that although PANDA can optionally 
take protein-protein interaction (PPI) information as an 
input, to generate a more fair comparison with other net¬ 
work reconstruction approaches PPI data was never used 
in any of the networks generated as a part of this study. 
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4.10. Calculating gene degree and comparing with 
DNase hypersensitivity data 


4.11. Clustering networks/expression and running 
GSEA 


When comparing with the DNase information, we cal¬ 
culated gene degree three different ways: 




C 1 = 

i 

4 m+e) = 5>W } 


(6) 


where is the “probability” that TF i is expressed, 
calculated by taking the inverse CDF of the z-score of 
the TF’s expression compared to the background of its 
expression in all other samples; is the “probability” 
of the edge from the LIONESS networks, found by tak¬ 
ing the inverse CDF of the predicted edge-weight score 
(which is in z-score units); is the “probability” of the 
edge from the motif data, either 0 or 1 based on whether 
the motif of TF i was found in the promoter of gene j. 


Clusters of networks and expression samples were 
defined based on a hierarchical clustering in which we 
row-normalized edge-weights (or gene expression) across 
samples, calculated distance based on the Spearman 
correlation, and performed a complete-linkage cluster¬ 
ing. When performing Gene Set Enrichment Analysis 
(GSEA) on the network targeting we calculated a score 
for each gene by summing up all predicted edge-weight 
scores for each gene. We then used LIMMA analysis m 
to determine the differences between the clusters based 
on these scores and ranked genes based on the log-fold 
changes. We ran a pre-ranked GSEA with 1000 iter¬ 
ations and observed 68 significant Reactome pathway 
signatures [32] with an FDR <0.1 and a gene set size less 
than forty. An equivalent clustering/LIMMA/GSEA 
analysis on the gene expression data resulted in no 
significant Reactome pathways. 
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5. SUPPLEMENTAL FIGURES 


Create 10000 networks: 

-> 40 nodes each 

-> same degree-distribution 

-> all nodes have both in/out degree 


Run Boolean model: 

-> 1000 random initial states 
-> 1000 corresponding steady states 
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Average steady state values to 

get expected "expression" levels 

for genes in each network 


Apply LIONESS: 

-> use Pearson correlation to 
estimate the aggregate 
network models 
-> Obtain 10000 single-sample 
networks 
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LIONESS-predicted Single Sample Networks 


Supplemental Figure 1: A schematic overview of how we generated in silico expression data for a set of known underlying gene 
regulatory network models. 
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Supplemental Figure 2: Results when applying LIONESS to yeast cell cycle data using various other aggregate reconstruction 
approaches. Although LIONESS is a generalizable approach, these results highlight the importance of selecting a robust 
underlying aggregate reconstruction algorithm when applying LIONESS. 
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Supplemental Figure 3: A hierarchical clustering of the 153 single-sample networks predicted for LIONESS in lymphoblastoid 
cell-lines. Close inspection reveals that single-sample networks derived from separate expression samples (technical replicates), 
but measuring the same cell-line, almost always cluster together. The AUC of these networks, when using ChIP-seq data for 
the GM12878 cell-line as a gold-standard, is shown directly below the clustering. 
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6. SUPPLEMENTAL MATERIAL 

6.A. Derivation to Find Regulatory Networks for Individual Samples in a Collection 

To begin, we assume that the value of a given edge (e^) from a transcription factor (i) to a gene ( j ) predicted 
by a network reconstruction algorithm using a collection of samples (a) is the linear combination of the value of that 
edge across networks specific to each of the input samples (e^), where represents the relative contribution of 
sample ( s ): 


N 


N 


4 a) = where £™7 = 1 


(7) 


S = 1 


S = 1 


Given this assumption, we can calculate two “aggregate” networks, one using all samples (e^), as described above, 
and the other using all but one of the samples 


N 


N 


o (cx-q) _ 


= where ^ = 1 


( 8 ) 


s/g s/g 

Now, subtracting these two “aggregate” network estimates we get: 


N 


e ij e ij 


= £< 
s= 1 


N 

c ij 

(9) 

s/g 



N 

N 


+ £»F4' > 

— £ W s a ~ q ^ e if 

(10) 

s#g 

s^q 


N 



+ I>3 a) -' 

4 a ~ q) )e% 

(11) 

s/g 




We can then solve for the network specific to a the single sample, : 


M = 

(a) 

Wn 


W t 


(a) 


N 

4 a) - 4r 9) + E(^ a_9) - «'i a) )ei; ) 

s/g 
N 

4 “’-£ 




s#g 


( 12 ) 

(13) 


If we then stipulate that the weights used to estimate eff are related to the weights used to estimate e-J ^ by a 
constant, = Cwi a ~ q \ which, given Equation^ implies, (7 = 1 — Wq a \ we can calculate the values of the edges 
in the single-sample network, , in terms of the two “aggregate” networks: 


M = 

ij (a) 


s^q 

4“ l -(i-»5“ , )4r ) 


N 




w { 


(a) 


, 0 ) 


J a ) _ 

e ij e ij 


+e, 


(a-q) 

ij 


(14) 

(15) 

(16) 


Finally, we can simplify this equation by (optionally) assuming that each sample is given equal weight (Wq 0 > = jr ): 

(17) 


= N 


„(«) _ „(«-<?) 
ij ij 


+e (a ~ q) 
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6.B. Application to Pearson Correlation 

To begin, we remember that the Pearson correlation (r) between two variables, X and Y can be defined as: 

1 N 

r = —-— W 

N - iZj 


Xj — X\ Yj —Y 


Sx 


Sy 


, where X = 1 and , EE Y J {X i - X) 


(18) 


We can then use the Pearson to calculate two “aggregate” correlation networks, one using all samples (resulting in 
r), and the other using all but one of the samples (resulting in r'): 


r' = —-— V' 
N — 9 


N 


X — X r \ Yj — Y r 


N - 2 


i^q 


S'x 




1 N 

|, where X' = ^Xj and S'> 

i^q 


\ 


^ f^X, - X'f 


N - 2 




Now, using the “LIONESS” equation for deriving single samples we see that: 


4f = W(e£ } - 4+ e Yj ) 


4y = N(r xy -r' xy ) + r' x 


= N 


xy) i »xy 
N 


N 




Xi-X\ Yi-Y 


Sx 


Sy 


~^—Y 

N - 9 


N — 2 


i^q 


Xi-X'\ Yj-Y 


i^q 


1 y ^ t Xi-X' \ I Yi-Y' 

+ Af _ 9 


N 

= N — 1 ^ 


Xi-X\ Yi-Y 


N — 1 Y^ 

~ N - 9 2^ 


N - 2 


i^q 


Xi-X'\ Yi-Y 


S' Y 


Qf 


N X 0 -X\ Y a -Y 


N- 1 \ Sx 


Sy 


N X (Xi-x\ (Yi-Y\ N-lX I Xi-X'\ I Yi-Y' 


N 


XL ^ 




Sx 


Sy 


iv — 1 \ 

_ N - 2 ^ 




5V 




JV X„-X\ Ya-Y 


N - 1 \ 5 


+ E 

- 


A 


A- 1 


Xi-X\ Yi-Y 


N-l\ I Xi-X'\ ( Yi-Y 
N -2 


We also note that (from equations [l8| and [l9| above): 


x= X x - 

i i * 


* 


X 

i^q 

N- 1 

X 




And consequently, 


= *' + -(*,-*') 


Xi-X = Xi-x' + L(x' - X ? ) 


(19) 

( 20 ) 
( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 


Thus when the difference between X g and X is much less than the total number of samples being considered (X), 
X' X X and S' x X Sx • This is most likely for large values of N in which case the above equation can be simplified 
to: 


r (<?) 

xy 


N X a -X\ Y a -Y 


N - 1 


S'x 


Sy 


(31) 


We can now see that, for Pearson-correlation derived “single-sample” networks, the linearity assumption holds, given 
the above limits on N and the relationship between X q and X'. Specifically, if we use equation pi] and average over 
all “single-sample” Pearson-derived networks, we see: 

N . N 


iy r « - 1 v 


Xj — X \ [Yj — Y 


Sx 


Sy 


(32) 


which is equal to r (equation [l8| . 



























































